#include "mloop.h"
#include "grid.h"
#include <ctype.h>

char resname[26][4] = {"ALA", "B", "CYS", "ASP", "GLU", "PHE", "GLY", "HIS", "ILE", "J", "LYS", "LEU", "MET", "ASN", "O", "PRO", "GLN", "ARG", "SER", "THR", "U", "VAL", "TRP", "X", "TYR", "Z"};

int restype[26] = {0, -1, 1, 2, 3, 4, 5, 6, 7, -1, 8, 9, 10, 11, -1, 12, 13, 14, 15, 16, -1, 17, 18, -1, 19, -1};

int judge_type (char name[4])
{
	int i, type;

	for (i = 0; i < 26; i++)
	{
		if (!strncmp (name, resname[i], 3))
		{
			type = restype[i];
			break;
		}
	}
	if (i == 26) /* Treat Unkown Residue as Alanine */
		type = 0;

	return type;
}

void read_rotamer()
{
	FILE *fp;
	char label[256], buf[256];
	int i, j, k, n, m, p;

	fp = fopen ("ROTAMER", "r");
	if (fp == NULL)
	{
		printf ("open ROTAMER failed\n");
		exit (-1);
	}
	do
	{
		fscanf (fp, "%s", label);
		if (!strcmp (label, "RESI"))
		{
			int a, b, c;
			char tmp[4];

			fscanf (fp, "%s%d%d%d", tmp, &a, &b, &c);
			for (i = 0; i < 26; i++)
			{
				if (!strcmp (tmp, resname[i]))
				{
					n = restype[i];
					break;
				}
			}
			sidechain[n].atom_num = a;
			sidechain[n].angle_num = b;
			sidechain[n].rotamer_num = c;
			sidechain[n].rotamer = calloc (sizeof (struct rotamer), sidechain[n].rotamer_num);
			p = -1;
		}
		else if (!strcmp (label, "ROTAMER"))
		{
			p++;
			for (i = 0; i < sidechain[n].angle_num; i++)
				fscanf (fp, "%f", &sidechain[n].rotamer[p].angle[i]);
			m = 0;
		}
		else if (!strcmp (label, "ATOM"))
		{
			fgets (buf, 256, fp);
			strncpy (sidechain[n].atom[m].atomname, buf + 9, 3);
			strncpy (sidechain[n].atom[m].resname, buf + 13, 3);
			sscanf (buf + 26, "%f%f%f", &sidechain[n].rotamer[p].coor[m][0], &sidechain[n].rotamer[p].coor[m][1], &sidechain[n].rotamer[p].coor[m][2]);
			m++;
		}
	}
	while (!feof (fp));
	fclose (fp);

	for (i = 0; i < 20; i++)
	{
		for (j = 0; j < sidechain[i].rotamer_num; j++)
		{
			for (k = 0; k < sidechain[i].atom_num; k++)
			{
				sidechain[i].rotamer[j].deal[k][0] = sidechain[i].rotamer[j].coor[k][0] - sidechain[i].rotamer[j].coor[1][0];
				sidechain[i].rotamer[j].deal[k][1] = sidechain[i].rotamer[j].coor[k][1] - sidechain[i].rotamer[j].coor[1][1];
				sidechain[i].rotamer[j].deal[k][2] = sidechain[i].rotamer[j].coor[k][2] - sidechain[i].rotamer[j].coor[1][2];
			}
			{
				float x12, y12, z12, x13, y13, z13;
				x12 = sidechain[i].rotamer[j].deal[0][0];
				y12 = sidechain[i].rotamer[j].deal[0][1];
				z12 = sidechain[i].rotamer[j].deal[0][2];
				x13 = sidechain[i].rotamer[j].deal[sidechain[i].atom_num - 1][0];
				y13 = sidechain[i].rotamer[j].deal[sidechain[i].atom_num - 1][1];
				z13 = sidechain[i].rotamer[j].deal[sidechain[i].atom_num - 1][2];

				sidechain[i].rotamer[j].deal[1][0] = y12 * z13 - y13 * z12;
				sidechain[i].rotamer[j].deal[1][1] = z12 * x13 - z13 * x12;
				sidechain[i].rotamer[j].deal[1][2] = x12 * y13 - x13 * y12;
			}
			for (k = 2; k < sidechain[i].atom_num - 1; k++)
			{
				float tmp[3];

				tmp[0] = sidechain[i].rotamer[j].deal[k][0] * sidechain[i].rotamer[j].deal[0][0] + sidechain[i].rotamer[j].deal[k][1] * sidechain[i].rotamer[j].deal[0][1] + sidechain[i].rotamer[j].deal[k][2] * sidechain[i].rotamer[j].deal[0][2];
				tmp[1] = sidechain[i].rotamer[j].deal[k][0] * sidechain[i].rotamer[j].deal[sidechain[i].atom_num - 1][0] + sidechain[i].rotamer[j].deal[k][1] * sidechain[i].rotamer[j].deal[sidechain[i].atom_num - 1][1] + sidechain[i].rotamer[j].deal[k][2] * sidechain[i].rotamer[j].deal[sidechain[i].atom_num - 1][2];
				tmp[2] = sidechain[i].rotamer[j].deal[k][0] * sidechain[i].rotamer[j].deal[1][0] + sidechain[i].rotamer[j].deal[k][1] * sidechain[i].rotamer[j].deal[1][1] + sidechain[i].rotamer[j].deal[k][2] * sidechain[i].rotamer[j].deal[1][2];
				memcpy (sidechain[i].rotamer[j].deal[k], tmp, sizeof (tmp));
			}
		}
	}
}

void paste_sidechain(float N[3], float CA[3], float C[3], struct candidate *p, int begin, int resno, int type, int rotamer_no, int flag)
{
	int i;
	float a_data[3*3], b_data[3*3];

	if (type == 5) return; /* gly */
	if (flag)
	{
		for (i = 2; i < sidechain[type].atom_num - 1; i++)
		{
			strcpy (loop.loopatom[begin + i - 2].name.resname, sidechain[type].atom[i].resname);
			strcpy (loop.loopatom[begin + i - 2].name.atomname, sidechain[type].atom[i].atomname);
			loop.loopatom[begin + i - 2].radius = judge_radius (loop.loopatom[begin + i - 2].name);
			loop.loopatom[begin + i - 2].status = judge_status (loop.loopatom[begin + i - 2].name);
			loop.loopatom[begin + i - 2].resno = resno;
			loop.loopatom[begin + i - 2].build_flag = 1;
		}
	}

	a_data[0] = N[0] - CA[0];
	a_data[1] = N[1] - CA[1];
	a_data[2] = N[2] - CA[2];
	a_data[3] = C[0] - CA[0];
	a_data[4] = C[1] - CA[1];
	a_data[5] = C[2] - CA[2];
	a_data[6] = a_data[1] * a_data[5] - a_data[4] * a_data[2];
	a_data[7] = a_data[2] * a_data[3] - a_data[5] * a_data[0];
	a_data[8] = a_data[0] * a_data[4] - a_data[3] * a_data[1];

	matrix_inverse (a_data, b_data);

	for (i = 2; i < sidechain[type].atom_num - 1; i++)
	{
		p->coor[begin + i - 2][0] = CA[0] + b_data[0] * sidechain[type].rotamer[rotamer_no].deal[i][0] + b_data[1] * sidechain[type].rotamer[rotamer_no].deal[i][1] + b_data[2] * sidechain[type].rotamer[rotamer_no].deal[i][2];
		p->coor[begin + i - 2][1] = CA[1] + b_data[3] * sidechain[type].rotamer[rotamer_no].deal[i][0] + b_data[4] * sidechain[type].rotamer[rotamer_no].deal[i][1] + b_data[5] * sidechain[type].rotamer[rotamer_no].deal[i][2];
		p->coor[begin + i - 2][2] = CA[2] + b_data[6] * sidechain[type].rotamer[rotamer_no].deal[i][0] + b_data[7] * sidechain[type].rotamer[rotamer_no].deal[i][1] + b_data[8] * sidechain[type].rotamer[rotamer_no].deal[i][2];
	}
}

void build_torsion(float a[3], float b[3], float c[3], float d[3], float length, float angle, float dihedral)
{
	float x12, y12, z12, x13, y13, z13, x23, y23, z23;
	float a_data[3*3], b_data[3];
	float dist1, dist2, dist3;
	float cos1, sin1;

	x12 = b[0] - c[0];
	y12 = b[1] - c[1];
	z12 = b[2] - c[2];
	x13 = a[0] - b[0];
	y13 = a[1] - b[1];
	z13 = a[2] - b[2];
	x23 = a[0] - c[0];
	y23 = a[1] - c[1];
	z23 = a[2] - c[2];
	dist1 = sqrtf (x12 * x12 + y12 * y12 + z12 * z12);
	dist2 = sqrtf (x13 * x13 + y13 * y13 + z13 * z13);
	dist3 = sqrtf (x23 * x23 + y23 * y23 + z23 * z23);
	cos1 = (dist1 * dist1 + dist2 * dist2 - dist3 * dist3) / (2.0 * dist1 * dist2);
	sin1 = sqrtf (1.0 - cos1 * cos1);

	/* limited the point in the plane defined by a b c */
	a_data[0] = y12 * z13 - y13 * z12;
	a_data[1] = z12 * x13 - z13 * x12;
	a_data[2] = x12 * y13 - x13 * y12;
	b_data[0] = 0.0;

	/* limited the angle of b-c-d */
	a_data[3] = x12;
	a_data[4] = y12;
	a_data[5] = z12;
	b_data[1] = length * dist1 * cosf (angle / RADIAN);

	/* limited the angle of a-b and c-d */
	a_data[6] = x13;
	a_data[7] = y13;
	a_data[8] = z13;
	b_data[2] = -length * dist2 * (cosf (angle / RADIAN) * cos1 - sinf (angle / RADIAN) * sin1);

	linear_solve (a_data, b_data, d);
	d[0] += c[0];
	d[1] += c[1];
	d[2] += c[2];
	rotate_atom (b, c, d, dihedral);
}

/* four points are coplaner, the angle of CA-C-N is 117.5, the angle of O=C-N is 121 */
void build_N(float CA[3], float C[3], float O[3], float N[3])
{
	float x12, y12, z12, x13, y13, z13;
	float a_data[3*3], b_data[3];

	x12 = CA[0] - C[0];
	y12 = CA[1] - C[1];
	z12 = CA[2] - C[2];
	x13 = O[0] - C[0];
	y13 = O[1] - C[1];
	z13 = O[2] - C[2];

	/* limited the point in the plane defined by CA, C, O */
	a_data[0] = y12 * z13 - y13 * z12;
	a_data[1] = z12 * x13 - z13 * x12;
	a_data[2] = x12 * y13 - x13 * y12;
	b_data[0] = 0.0;

	/* limited the angle of CA-C-N */
	a_data[3] = x12;
	a_data[4] = y12;
	a_data[5] = z12;
	b_data[1] = 1.52 * 1.33 * cosf (117.5 / RADIAN);

	/* limited the angle of O=C-N */
	a_data[6] = x13;
	a_data[7] = y13;
	a_data[8] = z13;
	b_data[2] = 1.23 * 1.33 * cosf (121.0 / RADIAN);

	linear_solve (a_data, b_data, N);
	N[0] += C[0];
	N[1] += C[1];
	N[2] += C[2];
}

/* four points are coplaner, the angle of CA-C=O is 121.5, the angle of O=C-N is 121 */
void build_O(float CA[3], float C[3], float N[3], float O[3])
{
	float x12, y12, z12, x13, y13, z13;
	float a_data[3*3], b_data[3];

	x12 = CA[0] - C[0];
	y12 = CA[1] - C[1];
	z12 = CA[2] - C[2];
	x13 = N[0] - C[0];
	y13 = N[1] - C[1];
	z13 = N[2] - C[2];

	/* limited the point in the plane defined by CA, C, N */
	a_data[0] = y12 * z13 - y13 * z12;
	a_data[1] = z12 * x13 - z13 * x12;
	a_data[2] = x12 * y13 - x13 * y12;
	b_data[0] = 0.0;

	/* limited the angle of CA-C=O */
	a_data[3] = x12;
	a_data[4] = y12;
	a_data[5] = z12;
	b_data[1] = 1.52 * 1.23 * cosf (121.5 / RADIAN);

	/* limited the angle of O=C-N */
	a_data[6] = x13;
	a_data[7] = y13;
	a_data[8] = z13;
	b_data[2] = 1.23 * 1.33 * cosf (121.0 / RADIAN);

	linear_solve (a_data, b_data, O);
	O[0] += C[0];
	O[1] += C[1];
	O[2] += C[2];
}

void build_CA(float CA1[3], float C[3], float N[3], float CA2[3], int type)
{
	if (type == 1)
		build_torsion (CA1, C, N, CA2, 1.45, 120.0, 0.0);
	else
		build_torsion (CA1, C, N, CA2, 1.45, 120.0, 180.0);
}

void build_loop(char *pdbfile, int segid, char *loopseq, int start_resno, int begin, int end)
{
	FILE *fp;
	char buf[256];
	int i, type;
	float last_CA[3], last_C[3], last_O[3], last_N[3], grid_tmp[3];

	fp = fopen (pdbfile, "r");
	while (fgets (buf, 256, fp))
	{
		int resno;
		char tmp[5];
		struct atomname atom;

		if (strncmp (buf, "ATOM", 4) != 0)
			continue;
		strncpy (tmp, buf + 22, 4);
		resno = atoi (tmp);

		if (resno == begin - 1)
		{
			sscanf (buf + 30, "%f%f%f", &current[0].coor[loop.loopnum][0], &current[0].coor[loop.loopnum][1], &current[0].coor[loop.loopnum][2]);
			if (!strncmp (buf + 13, "N  ", 3))
				memcpy (grid_tmp, current[0].coor[loop.loopnum], sizeof (float) * 3);
			else if (!strncmp (buf + 13, "CA ", 3))
				memcpy (last_CA, current[0].coor[loop.loopnum], sizeof (float) * 3);
			else if (!strncmp (buf + 13, "C  ", 3))
				memcpy (last_C, current[0].coor[loop.loopnum], sizeof (float) * 3);
			else if (!strncmp (buf + 13, "O  ", 3))
				memcpy (last_O, current[0].coor[loop.loopnum], sizeof (float) * 3);
			strncpy (loop.loopatom[loop.loopnum].name.atomname, buf + 13, 3);
			strncpy (loop.loopatom[loop.loopnum].name.resname, buf + 17, 3);
			loop.loopatom[loop.loopnum].type = judge_type (loop.loopatom[loop.loopnum].name.resname);
			loop.loopatom[loop.loopnum].radius = judge_radius (loop.loopatom[loop.loopnum].name);
			loop.loopatom[loop.loopnum].status = judge_status (loop.loopatom[loop.loopnum].name);
			loop.loopatom[loop.loopnum].resno = resno;
			loop.loopatom[loop.loopnum].build_flag = 0;
			loop.loopnum++;
		}
		else if (resno > begin - 1)
			break;
	}
	fclose (fp);

	for (i = begin; i <= end; i++)
	{
		float phi, psi;
		type = 0;
		if (islower(loopseq[i - begin]))
		{
			type = 1;
			loopseq[i - begin] = toupper(loopseq[i - begin]);
		}
		if (!isupper (loopseq[i - begin]) || restype[loopseq[i - begin] - 'A'] == -1)
		{
			printf ("Unknown residue: %c\n", loopseq[i - begin]);
			exit (-1);
		}
		loop.type[start_resno + i - begin] = restype[loopseq[i - begin] - 'A'];
		loop.start[start_resno + i - begin] = loop.loopnum;
		if (loop.loopnum + sidechain[loop.type[start_resno + i - begin]].atom_num + 1 >= LOOPNUM)
		{
			printf ("LOOPNUM is too small\n");
			exit (-1);
		}
		/* build N atom */
		build_N (last_CA, last_C, last_O, last_N);
		memcpy (current[0].coor[loop.loopnum], last_N, sizeof (float) * 3);
		strncpy (loop.loopatom[loop.loopnum].name.atomname, "N  ", 3);
		strncpy (loop.loopatom[loop.loopnum].name.resname, resname[loopseq[i - begin] - 'A'], 3);
		loop.loopatom[loop.loopnum].radius = judge_radius (loop.loopatom[loop.loopnum].name);
		loop.loopatom[loop.loopnum].status = judge_status (loop.loopatom[loop.loopnum].name);
		loop.loopatom[loop.loopnum].resno = i;
		loop.loopatom[loop.loopnum].build_flag = 1;
		loop.loopnum++;
		/* build CA atom */
		build_CA (last_CA, last_C, last_N, current[0].coor[loop.loopnum], type);
		memcpy (last_CA, current[0].coor[loop.loopnum], sizeof (float) * 3);
		strncpy (loop.loopatom[loop.loopnum].name.atomname, "CA ", 3);
		strncpy (loop.loopatom[loop.loopnum].name.resname, resname[loopseq[i - begin] - 'A'], 3);
		loop.loopatom[loop.loopnum].radius = judge_radius (loop.loopatom[loop.loopnum].name);
		loop.loopatom[loop.loopnum].status = judge_status (loop.loopatom[loop.loopnum].name);
		loop.loopatom[loop.loopnum].resno = i;
		loop.loopatom[loop.loopnum].build_flag = 1;
		loop.loopnum++;
		/* randomize mainchain dihedral (phi, psi) */
		rand_phi_psi (&phi, &psi, loop.type[start_resno + i - begin]);
		current[0].phi[start_resno + i - begin] = phi;
		current[0].psi[start_resno + i - begin] = psi;
		/* build C atom */
		build_torsion (last_C, last_N, last_CA, current[0].coor[loop.loopnum], 1.52, 111.6, phi);
		memcpy (last_C, current[0].coor[loop.loopnum], sizeof (float) * 3);
		strncpy (loop.loopatom[loop.loopnum].name.atomname, "C  ", 3);
		strncpy (loop.loopatom[loop.loopnum].name.resname, resname[loopseq[i - begin] - 'A'], 3);
		loop.loopatom[loop.loopnum].radius = judge_radius (loop.loopatom[loop.loopnum].name);
		loop.loopatom[loop.loopnum].status = judge_status (loop.loopatom[loop.loopnum].name);
		loop.loopatom[loop.loopnum].resno = i;
		loop.loopatom[loop.loopnum].build_flag = 1;
		loop.loopnum++;
		/* build O atom */
		build_torsion (last_N, last_CA, last_C, current[0].coor[loop.loopnum], 1.33, 117.5, psi);
		build_O (last_CA, last_C, current[0].coor[loop.loopnum], last_O);
		memcpy (current[0].coor[loop.loopnum], last_O, sizeof (float) * 3);
		strncpy (loop.loopatom[loop.loopnum].name.atomname, "O  ", 3);
		strncpy (loop.loopatom[loop.loopnum].name.resname, resname[loopseq[i - begin] - 'A'], 3);
		loop.loopatom[loop.loopnum].radius = judge_radius (loop.loopatom[loop.loopnum].name);
		loop.loopatom[loop.loopnum].status = judge_status (loop.loopatom[loop.loopnum].name);
		loop.loopatom[loop.loopnum].resno = i;
		loop.loopatom[loop.loopnum].build_flag = 1;
		loop.loopnum++;
		/* build sidechain */
		current[0].rotamer[start_resno + i - begin] = sidechain[loop.type[start_resno + i - begin]].rotamer_num * drand48();
		paste_sidechain (last_N, last_CA, last_C, &current[0], loop.loopnum, i, loop.type[start_resno + i - begin], current[0].rotamer[start_resno + i - begin], 1);
		loop.loopnum += sidechain[loop.type[start_resno + i - begin]].atom_num - 3;

		loop.end[start_resno + i - begin] = loop.loopnum - 1;
	}

	/* build dummy_N atom */
	build_N (last_CA, last_C, last_O, current[0].dummy_anchor[segid].N);
	/* build dummp_CA atom */
	if (islower(loopseq[i - begin]))
	{
		printf ("Anchor of %d loop is Proline.\n", segid);
		build_CA (last_CA, last_C, current[0].dummy_anchor[segid].N, current[0].dummy_anchor[segid].CA, 1);
	}
	else
		build_CA (last_CA, last_C, current[0].dummy_anchor[segid].N, current[0].dummy_anchor[segid].CA, 0);

	fp = fopen (pdbfile, "r");
	while (fgets (buf, 256, fp))
	{
		int resno;
		char tmp[5];
		struct atomname atom;

		if (strncmp (buf, "ATOM", 4) != 0)
			continue;
		strncpy (tmp, buf + 22, 4);
		resno = atoi (tmp);

		if (resno == end + 1)
		{
			if (loop.loopnum >= LOOPNUM)
			{
				printf ("LOOPNUM is too small\n");
				exit (-1);
			}
			sscanf (buf + 30, "%f%f%f", &current[0].coor[loop.loopnum][0], &current[0].coor[loop.loopnum][1], &current[0].coor[loop.loopnum][2]);
			if (!strncmp (buf + 13, "N  ", 3))
				memcpy (loop.anchor[segid].N, current[0].coor[loop.loopnum], sizeof (float) * 3);
			else if (!strncmp (buf + 13, "CA ", 3))
				memcpy (loop.anchor[segid].CA, current[0].coor[loop.loopnum], sizeof (float) * 3);
			else if (!strncmp (buf + 13, "C  ", 3))
			{
				grid_tmp[0] += current[0].coor[loop.loopnum][0];
				grid_tmp[1] += current[0].coor[loop.loopnum][1];
				grid_tmp[2] += current[0].coor[loop.loopnum][2];
			}
			strncpy (loop.loopatom[loop.loopnum].name.atomname, buf + 13, 3);
			strncpy (loop.loopatom[loop.loopnum].name.resname, buf + 17, 3);
			loop.loopatom[loop.loopnum].type = judge_type (loop.loopatom[loop.loopnum].name.resname);
			loop.loopatom[loop.loopnum].radius = judge_radius (loop.loopatom[loop.loopnum].name);
			loop.loopatom[loop.loopnum].status = judge_status (loop.loopatom[loop.loopnum].name);
			loop.loopatom[loop.loopnum].resno = resno;
			loop.loopatom[loop.loopnum].build_flag = 0;
			loop.loopnum++;
		}
		else if (resno > end + 1)
			break;
	}
	fclose (fp);

	grid.center[0] += grid_tmp[0] * 0.5;
	grid.center[1] += grid_tmp[1] * 0.5;
	grid.center[2] += grid_tmp[2] * 0.5;

#ifdef DEBUG
	//memcpy (&ref, &current[0], sizeof (struct candidate));
	fp = fopen (pdbfile, "r");
	while (fgets (buf, 256, fp))
	{
		int resno, m;
		char tmp[5];
		struct atomname atom;

		if (strncmp (buf, "ATOM", 4) != 0)
			continue;
		strncpy (tmp, buf + 22, 4);
		resno = atoi (tmp);

		if (resno >= begin && resno <= end)
		{
			if (!strncmp (buf + 13, "N  ", 3))
				m = loop.start[start_resno + resno - begin];
			else if (!strncmp (buf + 13, "CA ", 3))
				m = loop.start[start_resno + resno - begin] + 1;
			else if (!strncmp (buf + 13, "C  ", 3))
				m = loop.start[start_resno + resno - begin] + 2;
			else if (!strncmp (buf + 13, "O  ", 3))
				m = loop.start[start_resno + resno - begin] + 3;
			else if (!strncmp (buf + 13, "CB ", 3))
				m = loop.start[start_resno + resno - begin] + 4;
			else
				continue;

			sscanf (buf + 30, "%f%f%f", &(ref.coor[m][0]), &(ref.coor[m][1]), &(ref.coor[m][2]));
		}
	}
	fclose (fp);
#endif
}

void build_loops(char *pdbfile, int nseg, char loopseq[MAX_SEGMENT][20], int *begin, int *end)
{
	int i, j, k;
	double tmp1, tmp2; 

	read_rotamer ();

	loop.nseg = nseg;
	loop.loopnum = loop.resnum = loop.resnum2 = 0;
	grid.center[0] = grid.center[1] = grid.center[2] = 0.0;
	for (i = 0; i < nseg; i++)
	{
		build_loop (pdbfile, i, loopseq[i], loop.resnum, begin[i], end[i]);
		loop.seg_start[i] = loop.resnum;
		loop.resnum += end[i] - begin[i] + 1;
		loop.seg_end[i] = loop.resnum - 1;
	}
	grid.center[0] /= nseg;
	grid.center[1] /= nseg;
	grid.center[2] /= nseg;

	for (i = 0; i < loop.loopnum; i++)
	{
		if (loop.loopatom[i].build_flag == 0)
			continue;
	}

	for (i = 0; i < loop.resnum; i++)
	{
		for (j = loop.start[i]; j <= loop.end[i]; j++);
		{
			loop.loopatom[j].type = loop.type[i];
		}
	}

	for (i = 0; i < loop.loopnum; i++)
	{
		loop.loopatom[i].scale_factor = 1.0;
	}

	for (i = 0; i < loop.loopnum; i++)
	{
		if (loop.loopatom[i].build_flag == 0)
			continue;
		if (loop.loopatom[i].status & IS_BACKBONE)
			loop.loopatom[i].scale_factor = 1.0;
		else
			loop.loopatom[i].scale_factor = 0.75;
	}

	if (mode == 1)
		loop.probility = 2.0;
	else
	{
		tmp1 = tmp2 = 1.0;
		if (mode == 0)
		{
			for (i = 0; i < loop.resnum; i++)
			{
				if (sidechain[loop.type[i]].rotamer_num > 1)
				{
					loop.resno[loop.resnum2] = i;
					loop.resnum2 ++;
					tmp2 *= sidechain[loop.type[i]].rotamer_num;
				}
			}
		}
		else if (mode == 2)
		{
			for (i = 0; i < loop.nseg; i++)
			{
				for (j = loop.seg_start[i]; j <= loop.seg_end[i]; j++)
				{
					if (j - loop.seg_start[i] >= 2 && loop.seg_end[i] - j >= 2)
						continue;
					if (sidechain[loop.type[j]].rotamer_num > 1)
					{
						loop.resno[loop.resnum2] = j;
						loop.resnum2 ++;
						tmp2 *= sidechain[loop.type[j]].rotamer_num;
						for (k = loop.start[j]; k <= loop.end[j]; k++)
						{
							loop.loopatom[k].status |= IS_LOCAL_AA;
						}
					}
				}
			}
		}

		for (i = 0; i < loop.resnum; i++)
		{
			tmp1 *= 3.0;
		}

		if (loop.resnum2 == 0)
			loop.probility = 2.0;
		else
		{
			loop.probility = tmp1 / (tmp1 + tmp2);
			if (loop.probility < 0.1)
				loop.probility = 0.1;
			else if (loop.probility > 0.5)
				loop.probility = 0.5;
		}
	}
	printf ("Probility: %f\n", loop.probility);
}

void read_environatom (char *pdbfile, int nseg, int *begin, int *end)
{
	FILE *fp;
	char buf[256];
	float maxx, maxy, maxz;
	float minx, miny, minz;
	int i, flag;

	maxx = grid.center[0] + GRID_SIZE * GRIDMAX_X * 0.5 + 2 * R_MAX;
	maxy = grid.center[1] + GRID_SIZE * GRIDMAX_Y * 0.5 + 2 * R_MAX;
	maxz = grid.center[2] + GRID_SIZE * GRIDMAX_Z * 0.5 + 2 * R_MAX;
	minx = grid.center[0] - GRID_SIZE * GRIDMAX_X * 0.5 - 2 * R_MAX;
	miny = grid.center[1] - GRID_SIZE * GRIDMAX_Y * 0.5 - 2 * R_MAX;
	minz = grid.center[2] - GRID_SIZE * GRIDMAX_Z * 0.5 - 2 * R_MAX;

	fp = fopen (pdbfile, "r");
	if (fp == NULL)
	{
		printf ("open %s failed\n", pdbfile);
		exit (-1);
	}
	atomnum = 0;
	while (fgets (buf, 256, fp))
	{
		int resno;
		char tmp[5];
		float x, y, z;

		if (strncmp (buf, "ATOM", 4) != 0)
			continue;
		strncpy (tmp, buf + 22, 4);
		resno = atoi (tmp);
		flag = 1;
		for (i = 0; i < nseg; i++)
		{
			if (resno >= begin[i] - 1 && resno <= end[i] + 1)
			{
				flag = 0;
				break;
			}
		}
		if (flag == 1)
		{
			sscanf (buf + 30, "%f%f%f", &x, &y, &z);
			if (x > maxx || x < minx || y > maxy || y < miny || z > maxz || z < minz)
				continue;
			if (atomnum == ATOMNUM)
			{
				printf ("ATOMNUM is too small\n");
				exit (-1);
			}
			environatom[atomnum].resno = resno;
			environatom[atomnum].coor[0] = x;
			environatom[atomnum].coor[1] = y;
			environatom[atomnum].coor[2] = z;
			strncpy (environatom[atomnum].name.atomname, buf + 13, 3);
			strncpy (environatom[atomnum].name.resname, buf + 17, 3);
			environatom[atomnum].type = judge_type (environatom[atomnum].name.resname);
			environatom[atomnum].radius = judge_radius (environatom[atomnum].name);
			environatom[atomnum].status = judge_status (environatom[atomnum].name);
			atomnum++;
		}
	}
	fclose (fp);
}
